##################
####One by One####
##################

#Figure A.6

#Plots for #Model 33 to #Model 39

##########
##########
#Model 33#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 1)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", 
                        
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo1.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist",
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 34#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 2)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil",
                        
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo2.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 35#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 3)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil","lag.lngas", 
                        
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil+lag.lngas
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo3.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 36#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 4)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil","lag.lngas", "lag.lnbdist1",
                         
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil+lag.lngas+lag.lnbdist1
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo4.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 37#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) +
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 5)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil","lag.lngas", "lag.lnbdist1",
                        "lag.lnpop_gpw_sum", 
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil+lag.lngas+lag.lnbdist1+lag.lnpop_gpw_sum
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo5.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 38#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 6)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil","lag.lngas", "lag.lnbdist1",
                        "lag.lnpop_gpw_sum","lag.lnexcluded",  
                        "ccode", "year", "gid"
)])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil+lag.lngas+lag.lnbdist1+lag.lnpop_gpw_sum+lag.lnexcluded
               
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo6.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

##########
##########
#Model 39#
##########
##########
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + 
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (Model 7)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c( "lnNTL","lag.cost_sanction","lag.lncapdist",
                        "lag.p_polity2", "lag.lnoil","lag.lngas", "lag.lnbdist1",
                        "lag.lnpop_gpw_sum","lag.lnexcluded", "lag.lngcp", 
                        "ccode", "year", "gid"
                        )])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.cost_sanction+lag.lncapdist+lag.cost_sanction:lag.lncapdist+
                 lag.p_polity2+lag.lnoil+lag.lngas+lag.lnbdist1+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lngcp
                 
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("obo7.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.cost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

